/** @file pra_matrix.c
@brief matrix opration functions
this file cantains all matrix oprations related to matrix
@author Prabhakar Pujeri
@bug No know bugs.
 */

#include "pramatrix.h"
#include <stdio.h>

pra_matrix* pra_matrix_alloc ( int row, int col )
{
    pra_matrix *t;
    t = ( pra_matrix * ) malloc ( sizeof ( pra_matrix ) );
    t->matrix = 0x0;
    t->col = col;
    t->row = row;
    t->matrix = ( double * ) malloc ( ( sizeof ( double ) *row*col ) );
    return t;
}

void pra_matrix_free ( pra_matrix *A)
{
    free ( A->matrix );
    free( A );
}

pra_matrix* pra_matrix_input ( pra_matrix *t )
{
    int i,j;
    if ( t->matrix > 0 )
    {
        for ( i = 0; i < t->row; i++ ) {
            for ( j = 0; j< t->col; j++ ) {
                printf ( "enter value for %dth row and %dth col\n",i,j );
                scanf ( "%lf",&t->matrix[(t->col*i)+j] );
            }
        }
    }
    return t;
}

void pra_matrix_output ( pra_matrix *t )
{
    int i,j;
    if ( t->matrix > 0 )
    {
        for ( i = 0; i < t->row; i++ ) {
            for ( j = 0; j< t->col; j++ ) {
                printf ( "%lf\t",(t->matrix[(t->col*i)+j]) );
            }
            printf ( "\n" );
        }
    }
}

void pra_matrix_add ( pra_matrix *A, pra_matrix *B, pra_matrix *C ) 
{
    int i,j;
    if ( (A->row == B->col) && (A->matrix > 0) && (B->matrix > 0) ) {
        for ( i = 0; i < A->row; i++ ) {
            for ( j = 0; j< B->col; j++ ) {
                C->matrix[(C->col*i)+j] = A->matrix[(A->col*i)+j]+B->matrix[(B->col*i)+j];
            }
        }
    } 
}
void pra_matrix_add_const ( pra_matrix *A, double n, pra_matrix *C )
{

}
void pra_matrix_sub ( pra_matrix *A, pra_matrix *B, pra_matrix *C ) 
{
    int i,j;
    if((A->row == B->col) && (A->matrix > 0) && (B->matrix > 0)) {
        for ( i = 0; i < A->row; i++ ) {
            for ( j = 0; j< B->col; j++ ) {
                C->matrix[(C->col*i)+j] = A->matrix[(A->col*i)+j]- B->matrix[(B->col*i)+j];
            }
        }
    } 
}
void pra_matrix_sub_const ( pra_matrix *A, double n, pra_matrix *C )
{
}
void pra_matrix_multi ( pra_matrix *A ,pra_matrix *B,pra_matrix *C )
{
    int i,j,k;
    if ((A->matrix > 0) && (B->matrix > 0)) {
        if ( B->col == A->row ) {
            for ( i = 0; i < A->row; i++ ) {
                for ( j = 0; j< B->col; j++ ) {
                    for( k = 0; k < A->col; k++ ) {
                        C->matrix[(A->col*i)+j] += A->matrix[(A->col*i)+k] * B->matrix[(k*A->col)+j];
                    }
                }
            }
        }
    }
}


void pra_matrix_dot_product ( pra_matrix *A, pra_matrix *B, pra_matrix *C ) 
{
    int i,j;
    if ( (A->matrix > 0) && (B->matrix > 0) ) {
        for ( i = 0; i < A->row; i++ ) {
            for ( j = 0; j< B->col; j++ ) {
                C->matrix[(C->col*i)+j] = A->matrix[(A->col*i)+j] * B->matrix[(B->col*i)+j];
            }
        }
    }
}

void pra_matrix_scalar_product ( pra_matrix *A, double n, pra_matrix *B ) 
{
    int i,j;
    if ( A->matrix > 0 ) {
        for ( i = 0; i < A->row; i++ ) {
            for ( j = 0; j< A->col; j++ ) {
                B->matrix[(A->col*i)+j] = n*(double)A->matrix[(A->col*i)+j] ;
            }
        }
    }
}
double pra_matrix_determinant ( pra_matrix *A, int n )
{
    int i,j,j1,j2;
    double det = 0;
    pra_matrix *M = NULL;

    if (n < 1) {

    } else if (n == 1) { 
        det = A->matrix[(0*n)+0];
    } else if (n == 2) {
        det = A->matrix[(0*n)+0] * A->matrix[(1*n)+1] - A->matrix[(1*n)+0] * A->matrix[(0*n)+1];
    } else {
        det = 0;
        for (j1 = 0;j1<n;j1++) {
            M = pra_matrix_alloc(n-1,n-1);
            for (i = 0;i<n-1;i++)
            for (i = 1;i<n;i++) {
                j2 = 0;
                for (j=0;j<n;j++) {
                    if (j == j1)
                    continue;
                    M->matrix[((i-1)*(n-1))+j2] = A->matrix[(i*n)+j];
                    j2++;
                }
            }
            det += pow(-1.0,j1+2.0) * A->matrix[(0*n)+j1] * pra_matrix_determinant(M,n-1);
            pra_matrix_free ( M );
        }
    }
    return(det);
}

void pra_matrix_transpose ( pra_matrix *A, int n, pra_matrix *T )
{
    int i,j;
    for (i=0;i<n;i++) {
        for (j=0;j<n;j++) {
            T->matrix[(i*n)+j] = A->matrix[(j*n)+i];
            T->matrix[(j*n)+i] = A->matrix[(i*n)+j];
        }
    }
}

void pra_matrix_co_factor ( pra_matrix *A, int n, pra_matrix *B )
{
    int i,j,ii,jj,i1,j1;
    double det;
    pra_matrix *C;
    C = pra_matrix_alloc(n-1,n-1);
    for (j=0;j<n;j++) {
        for (i=0;i<n;i++) {

            /* Form the adjoint a_ij */
            i1 = 0;
            for (ii=0;ii<n;ii++) {
                if (ii == i)
                continue;
                j1 = 0;
                for (jj=0;jj<n;jj++) {
                    if (jj == j)
                    continue;
                    C->matrix[(i1*(n-1))+j1] = A->matrix[(ii*n)+jj];
                    j1++;
                }
                i1++;
            }

            /* Calculate the determinate */
            det = pra_matrix_determinant (C,n-1);

            /* Fill in the elements of the cofactor */
            B->matrix[(i*n)+j] = pow(-1.0,i+j+2.0) * det;
        }
    }
    pra_matrix_free ( C );
}

void pra_matrix_adjoint ( pra_matrix *A, pra_matrix *B )
{
    pra_matrix *Co;
    Co=pra_matrix_alloc ( A->row, A->col );
    pra_matrix_co_factor ( A, A->row, Co );
    pra_matrix_transpose ( Co, A->row, B );
}

void pra_matrix_inverse ( pra_matrix *A, pra_matrix *B ) 
{
    double det;
    pra_matrix *T,*Co;
    Co = pra_matrix_alloc ( A->row, A->col );
    T = pra_matrix_alloc ( A->row, A->col );
    pra_matrix_co_factor ( A, A->row, Co );
    pra_matrix_transpose ( Co, A->row, T );
    det = pra_matrix_determinant ( A, A->row );
    pra_matrix_scalar_product ( T, (1/det), B );
}

int pra_matrix_equal ( pra_matrix *A, pra_matrix *B )
{
    int i,j,ret=-1;
    if((A->row == B->col) && (A->matrix > 0) && (B->matrix > 0)) {
        for ( i = 0; i < A->row; i++ ) {
            for ( j = 0; j< B->col; j++ ) {
                if ( A->matrix[(A->col*i)+j] == B->matrix[(B->col*i)+j] ) {
                    ret = 0;
                    continue;
                }else {  
                    ret = 1;
                    break;
                    break;
                }
            }
        }
    } 
    return ret;
}
int pra_matrix_symmetric( pra_matrix *A )
{
    int eq;
    pra_matrix *temp;
    temp = pra_matrix_alloc ( A->row, A->col );
    pra_matrix_transpose ( A, A->row, temp );
    eq=pra_matrix_equal ( A, temp );
    return eq;
}
pra_matrix* identity_matrix ( int n )
{
    pra_matrix *I;
    int i;
    I=pra_matrix_alloc(n,n);
    for ( i=0; i < n; i++ ) {
    }
    return I;
}

void  pra_matrix_set_to_null (pra_matrix *t)
{
    int i,j;
    if ( t->matrix > 0)
    {
        for ( i=0; i < t->row; i++ ) {
            for ( j=0; j< t->col; j++ ) {
                t->matrix[(t->col*i)+j]=0;
            }
            printf("\n");
        }
    }
}
